*Purpose: Post-Matching Estimation with Difference-in-Difference
** Main specifications for manuscript tables, including binary and continuous treatment ** 

clear


ssc install schemepack
set scheme white_tableau, perm

**0. Load data
use "$data\MASTER_MATCHED_R2.dta", clear
rename cell_uniqueID cellID

**1. Define globals 
global invarcontrols elevation slope NearCity_dist_km NearCity_population i.LandSurface_code i.Bioclimate_code i.Geology_code 
global timevarcontrols Precipitation Temperature 

** 2. Estimations
**DID with binary and continuous treatment 
*(1) DID 
reg NetOverNF treat post treat_post $invarcontrols $timevarcontrols , cluster(cd_mun)
estimates store col1
estadd local Mun_FE "No": col1
estadd local Cell_FE "No": col1
*(2) DID + Municipality FE
xtreg NetOverNF treat post treat_post $invarcontrols $timevarcontrols , fe i(cd_mun) cluster(cd_mun)
estimates store col2
estadd local Mun_FE "Yes": col2
estadd local Cell_FE "No": col2
*(3) DID + Cell FE
xtreg NetOverNF post treat_post $timevarcontrols, fe i(cellID) cluster(cd_mun)
estimates store col3
estadd local Mun_FE "No": col3
estadd local Cell_FE "Yes": col3

*(1) DID continuous
reg NetOverNF treat_percent post post_treat_percent $invarcontrols $timevarcontrols  , cluster(cd_mun)
estimates store col4
estadd local Mun_FE "No": col4
estadd local Cell_FE "No": col4
*(2) DID + Municipality FE continuous 
xtreg NetOverNF treat_percent post post_treat_percent $invarcontrols $timevarcontrols , fe i(cd_mun) cluster(cd_mun)
estimates store col5
estadd local Mun_FE "Yes": col5
estadd local Cell_FE "No": col5
*(3) DID + Cell FE continuous 
xtreg NetOverNF post post_treat_percent $timevarcontrols, fe i(cellID) cluster(cd_mun)
estimates store col6
estadd local Mun_FE "No": col6
estadd local Cell_FE "Yes": col6

*Print table in stata
esttab col1 col2 col3 col4 col5 col6 , keep(treat_post treat post_treat_percent treat_percent post) order(treat_post treat post_treat_percent treat_percent post) mtitles("Bin" "BinMunFE" "BunCellFE" "Cont" "ContMunFE" "ContCellFE")

*LaTeX Output Table
esttab  col1 col2 col3 col4 col5 col6  using "$tables\Table1.tex", label replace  starlevels(* 0.10 ** 0.05 *** 0.01)  se nonotes ///
keep(treat_post treat post_treat_percent treat_percent post) order(treat_post treat post_treat_percent treat_percent post) ///
varwidth(3)  style(tex) mtitles("" "" "" "" "" "" ) b(%12.3f) se(%12.3f) scalar("Mun_FE Mun FE" "Cell_FE Cell FE" )  ///
 prehead({ \begin{tabular}{lcccccc}    ///
	 \hline \hline ///
	 & \multicolumn{3}{c}{Binary treatment} & \multicolumn{3}{c}{Proportional treatment} \\)  ///


quietly: reg NetOverNF treat_percent post post_treat_percent `invarcontrols' `timevarcontrols'  , cluster(cd_mun)
estimates store ols, title(OLS)
quietly: xtreg NetOverNF treat_percent post post_treat_percent `invarcontrols' `timevarcontrols', fe i(cd_mun) cluster(cd_mun)
estimates store munfe, title(Mun FE)
quietly: xtreg NetOverNF post post_treat_percent `timevarcontrols', fe i(cellID) cluster(cd_mun)	
estimates store cellfe, title(Cell FE)
	
*Figure appendix with all FE.
coefplot (ols, transform(* = 0.25*(@)) label(OLS) ciopts(lwidth(0.75) lcolor("51 166 166")) mcolor("51 166 166") msymbol(Oh) msize(large)) (munfe, transform(* = 0.25*(@)) label(Mun FE) mcolor("242 149 68") msymbol(O) msize(large) ciopts(lwidth(0.75) lcolor("242 149 68")))  (cellfe, transform(* = 0.25*(@)) label(Cell FE) mcolor("242 80 65") msymbol(Sh) msize(large) ciopts(lwidth(0.75) lcolor("242 80 65")))  ///
 (ols, transform(* = 0.5*(@)) label("") mcolor("51 166 166") msymbol(Oh) msize(large) ciopts(lwidth(0.75) lcolor("51 166 166"))) (munfe, transform(* = 0.50*(@)) label("") mcolor("242 149 68") msymbol(O) msize(large) ciopts(lwidth(0.75) lcolor("242 149 68"))) ///
 (cellfe, transform(* = 0.50*(@)) label("") mcolor("242 80 65") msymbol(Sh) msize(large) ciopts(lwidth(0.75) lcolor("242 80 65")))  ///
(ols, transform(* = 0.75*(@)) label("")  mcolor("51 166 166") msymbol(Oh) msize(large) ciopts(lwidth(0.75) lcolor("51 166 166"))) (munfe, transform(* = 0.75*(@)) label("")  mcolor("242 149 68") msymbol(O) msize(large) ciopts(lwidth(0.75) lcolor("242 149 68"))) ///
 (cellfe, transform(* = 0.75*(@)) label("") mcolor("242 80 65") msymbol(Sh) msize(large) ciopts(lwidth(0.75) lcolor("242 80 65")))  ///
(ols, label("")  mcolor("51 166 166") msymbol(Oh) msize(large) ciopts(lwidth(0.75) lcolor("51 166 166"))) (munfe, label("")  mcolor("242 149 68") msymbol(O) msize(large) ciopts(lwidth(0.75) lcolor("242 149 68")))  (cellfe, label("") mcolor("242 80 65") msymbol(Sh) msize(large) ciopts(lwidth(0.75) lcolor(gs8))) , /// 
vertical yscale(range(0 0.3)) ylabel(0(.1)0.30) group(2 4 6 = "25%") keep(post_treat_percent) ytitle("Marginal effect on forest restoration")  coeflabels(post_treat_percent=" ") ///
 legend(pos(6) order(2 4 6) rows(1))  title("")  b1title("     25%                          50%                        75%                         100%") xtitle("") xline(0.77 1 1.23)
graph export "$figures\IntensityWFE.eps", as(eps) replace

*Figure with just cell fe

coefplot (cellfe, transform(* = 0.25*(@)) label("25%") mlabel(@plot) mcolor("242 80 65") msymbol(O) msize(large) ciopts(lwidth(0.75) lcolor("51 166 166")))  ///
 (cellfe, transform(* = 0.50*(@)) label("50%") mlabel(@plot)   mcolor("242 80 65") msymbol(O) msize(large) ciopts(lwidth(0.75) lcolor("51 166 166")))  ///
 (cellfe, transform(* = 0.75*(@)) label("75%") mlabel(@plot)  mcolor("242 80 65") msymbol(O) msize(large) ciopts(lwidth(0.75) lcolor("51 166 166")))  ///
 (cellfe, label("100%") mlabel(@plot)  mcolor("242 80 65") msymbol(O) msize(large) ciopts(lwidth(0.75) lcolor("51 166 166"))) , coeflabels(post_treat_percent=" ") /// 
vertical ylabel(0(.1)0.30, noticks nolabel) keep(post_treat_percent) ytitle(" ")   ///
 legend(pos(6) order(2 1) label(1 "95% CI") label(2 "Estimated effect by percentage of cell in restoration") rows(1))  title("")  ///
 xtitle("") saving(cellfe, replace) xlabel( , noticks nolabel)  plotregion(margin(zero)) 
 
*(3) DID + Cell FE
quietly: xtreg NetOverNF post treat_post `timevarcontrols' , fe i(cellID) cluster(cd_mun)
est store main
coefplot (main, label("Binary effect") ciopts(lwidth(0.75) lcolor("140 31 102")) msymbol(O) msize(large) mcolor("50 54 115")),  ///
	vertical  yscale(range(0 0.3)) ylabel(0(.1)0.30)  keep(treat_post) coeflabels(treat_post="")  xtitle("")  plotregion(margin(zero)) ///
	  ytitle("Marginal effect on restored forest/initial non-forest", size(medium)) legend( pos(6) order(2) label(1 "Binary effect")  rows(1))   ///
	xlabel( , noticks nolabel) scheme(white_tableau)   saving(basicTE, replace) fxsize(25) fysize(100)
	
gr combine basicTE.gph cellfe.gph, imargin(0 0 0 0)  ycommon col(2)
graph export "$figures\Figure4.eps", as(eps) replace


*Estimating total hectares restored

egen totalnf = sum(COV_NONFOREST) if treat == 1 & post == 0
sum totalnf
scalar baseline = r(mean)

gen mfx = 0.208*treat_percent*COV_NONFOREST if treat == 1 & post == 0
egen totalrecover = sum(mfx)
sum totalrecover
scalar recover = r(mean)
scalar list baseline recover

